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We introduce the heat method for computing the shortest geodesic distance 
to a specified subset {e.g., point or curve) of a given domain. The heat 
method is robust, efficient, and simple to implement since it is based on 
solving a pair of standard linear elliptic problems. The method represents a 
significant breakthrough in the practical computation of distance on a wide 
variety of geometric domains, since the resulting linear systems can be pref- 
actored once and subsequently solved in near-linear time. In practice, dis- 
tance can be updated via the heat method an order of magnitude faster than 
with state-of-the-art methods while maintaining a comparable level of accu- 
racy. We provide numerical evidence that the method converges to the exact 
geodesic distance in the limit of refinement; we also explore smoothed ap- 
proximations of distance suitable for applications where more regularity is 
required. 

Categories and Subject Descriptors: 1.3.7 [Computer Graphics]: Com- 
putational Geometry and Object Modeling — Geometric algorithms, lan- 
guages, and systems 

General Terms: digital geometry processing, discrete differential geometry, 
geodesic distance, distance transform, heat kernel 



1. INTRODUCTION 

Imagine touching a scorching hot needle to a single point on a sur- 
face. Over time heat spreads out over the rest of the domain and can 
be described by a function kt,x{y) called the heat kernel, which 
measures the heat transferred from a source x to a destination y af- 
ter time t. A well-known relationship between heat and distance is 
Varadhan's formula 1 1967], which says that the geodesic distance 
(j) between any pair of points x, 2/ on a Riemannian manifold can be 
recovered via a simple pointwise transformation of the heat kernel: 



= \im J -4:t log kt,x{y)- (1) 
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Fig. 1 . Geodesic distance on the Stanford Bunny. The heat method allows 
distance to be rapidly updated for new source points or curves. 

The intuition behind this behavior stems from the fact that heat dif- 
fusion can be modeled as a large collection of hot particles taking 
random walks starting at x: any particle that reaches a distant point 
y after a small time t has had little time to deviate from the short- 
est possible path. To date, however, this relationship has not been 
exploited by numerical algorithms that compute geodesic distance. 

Why has Varadhan's formula been overlooked in this context? 
The main reason, perhaps, is that it requires a precise numerical 
reconstruction of the heat kernel, which is difficult to obtain for 
small values of t - applying the formula to a mere approximation of 
kt,x does not yield the correct result, as illustrated in Figures[2]and 
[6] The main idea behind the heat method is to circumvent this issue 
by working with a broader class of inputs, namely any function 
whose gradient is parallel to geodesies. We can then separate the 
computation of distance into two separate stages: first compute the 
gradient of the distance field, then recover the distance itself. 

Relative to existing algorithms, the heat method offers two major 
advantages. First, it can be applied to virtually any type of geomet- 
ric discretization, including regular and irregular grids, polygonal 
meshes, and even unstructured point clouds. Second, it involves 
only the solution of sparse linear systems, which can be prefac- 
tored once and rapidly re- solved many times. This feature makes 
the heat method particularly valuable for applications such as shape 
matching, path planning, and level set-based simulation (e.g., free- 
surface fluid flows), which require repeated distance queries on a 
fixed geometric domain. Moreover, because linear elliptic equa- 
tions are widespread in scientific computing, the heat method can 
immediately take advantage of new developments in numerical lin- 
ear algebra and parallelization. 
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Fig. 2. Given an exact reconstruction of the heat kernel (top left) Varad- 
han's formula can be used to recover geodesic distance (bottom left) but 
fails in the presence of approximation or numerical error (middle, right), 
as shown here for a point source in ID. The robustness of the heat method 
stems from the fact that it depends only on the direction of the gradient. 

2. RELATED WORK 

The prevailing approach to distance computation is to solve the 
eikonal equation 

|V(/)| = 1 (2) 

subject to boundary conditions = over some subset 7 of 
the domain. This formulation is nonlinear and hyperbolic, mak- 
ing it difficult to solve directly. Typically one applies an iterative 
relaxation scheme such as Gauss-Seidel - special update orders 
are known as fast marching and fast sweeping, which are some 
of the most popular algorithms for distance computation on reg- 
ular grids ISethian 1996] and triangulated surfaces [Kimme l and| 
[Sethian 199 8 1. These algorithm s can also be used on im plicit sur- 
faces [Me moli and Sapiro 2001), point c louds | Memoli and Sapiro| 
^005], and polygon soup |Campen and Kobbelt 2011] , but only 
indirectly: distance is computed on a simplicial mesh or regular 
grid that approximates the original domain. Implementation of fast 
marching on simplicial grids is challenging due to the need for 
nonobtuse triangulations (which are notoriously difficult to obtain) 
or else a complex unfolding procedure to preserve monotonicity 
of the solution; moreover these issues are not well-studied in di- 
mensions greater than two. Fast marching and fast sweeping have 
asymptotic complexity of O(nlogn) and 0(n), respectively, but 
sweeping is often slower due to the lar ge number of sw eeps re- 
quired to obtain accurate results |Hysing and Turek 2005| . 

The main drawback of these methods is that they do not reuse 
information: the distance to different subsets 7 must be computed 
entirely from scratch each time. Also note that both sweeping and 
marching present challenges for parallelization: priority queues are 
inherently serial, and irregular meshes lack a natural sweeping or- 
der. Weber etal |2008 1 address this issue by decomposing surfaces 
into regular grids, but this decomposition resamples the surface and 
requires a low-distortion parameterization that may be difficult to 
obtain (note that the heat method would also benefit from such a 
decomposition). 

In a different development, Mitchell et al. 1| 1987| give an 
0(n^ log n) algorithm for computing the exact polyhedral distance 
from a single source to all other vertices of a triangulated surface. 
Surazhsky et al. | 2005 1 demonstrate that this algorithm tends to 
run in sub-quadratic time in practice, and present an approximate 
0{n log n) version of the algorithm with guaranteed error bounds; 
Bommes and Kobbelt |2007 | extend the algorithm to polygonal 
sources. Similar to fast marching, these algorithms propagate dis- 
tance information in wavefront order using a priority queue, again 
making them difficult to parallelize. More importantly, the amor- 




Fig. 3. The heat method computes the shortest distance to a subset 7 of a 
given domain. Gray curves indicate isolines of the distance function. 

tized cost of these algorithms (over many different source subsets 
7) is substantially greater than for the heat method since they do 
not reuse information from one subset to the next. Finally, although 
[Surazhsky et al. 2005) greatly simplifies the original formulation, 
these algorithms remain challenging to implement and do not im- 
mediately generalize to domains other than triangle meshes. 

Closest to our approach is the recent method of Rangarajan 
and Gurumoorthy |2011|, who do not appear to be aware of 
Varadahn's formula - they instead derive an analogous relationship 
(j) — —yjh log-^ between the distance function and solutions '0 to 
the time-independent Schrodinger equation. We emphasize, how- 
ever, that this derivation applies only in where takes a special 
form - in this case it may be just as easy to analytically invert the 
Euclidean heat kernel Ut,x = (47rt)~^/^e~'^*^'^'^^^/^*. Moreover, 
they compute solutions using the fast Fourier transform, which lim- 
its computation to regular grids. To obtain accurate results their 
method requires either the use of arbitrary-precision arithmetic or a 
combination of multiple solutions for various values of h\ no gen- 
eral guidance is provided for determining appropriate val ues of h. 

Finally, there is a large literature on smooth distance s |Coifman| 
[and Lafon 2'0Q6l|Fouss et al. 2007[[Lipman et al. 2010| , which are 
valuable in contexts where differentiability is required. However, 
existing smooth distances may not be appropriate in contexts where 
the geometry of the original domain is important, since they do not 
attempt to approximate the original metric and therefore substan- 
tially violate the unit-speed nature of geodesies (Figure[TO|. Inter- 
estingly enough, these distances also have an interpretation in terms 
of simple discretizations of heat flow - see Section [33] for further 
discussion. 




Fig. 4. Distance to the boundary on a region in the plane (left) or a surface 
in is achieved by simply placing heat along the boundary curve. Note 
good recovery of the cut locus, i.e., points with more than one closest point 
on the boundary. 
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Fig. 5. Outline of the heat method. (I) Heat u is allowed to diffuse for a 
brief period of time {left). (II) The temperature gradient Vw {center left) is 
normalized and negated to get a unit vector field X {center right) pointing 
along geodesies. (Ill) A function cf) whose gradient follows X recovers the 
final distance {right). 

3. THE HEAT METHOD 

Our method can be described purely in terms of operations on 
smooth manifol ds; we expl ore spatial and temporal discretiza- 
tion in Sections |3.1| and |3.2| respectively. Let A be the negative- 
semidefinite Laplace-Beltrami operator acting on (weakly) differ- 
entiable real-valued functions over a Riemannian manifold (M, g). 
The heat method consists of three basic steps: 

Algorithm 1 The Heat Method 

L Integrate the heat flow u — Au for some fixed time t. 
11. Evaluate the vector field X = — Vi^/| Vi^|. 
III. Solve the Poisson equation Acj) = V • X. 



The function approximates geodesic distance, approaching the 
true distance as t goes to zero (Eq. ([TJ). Note that the solution to 
step III is unique only up to an additive constant - final values sim- 
ply need to be shifted such that the smallest distance is zero. Initial 
conditions uq = 6{x) (i.e., a Dirac delta) recover the distance to a 
single source point x G M as in Figure [T] but in general we can 
compute the distance t o any pi ecewise submanifold 7 by setting uq 
to a generalized Dirac | | Villa 20 06] over 7 (see Figures[3]and|4]). 

The heat method can be motivated as follows. Consider an ap- 
proximation Ut of heat flow for a fixed time t. Unless Ut ex- 
hibits precisely the right rate of decay, Varadhan's transformation 
Ut ^ \/—4:t log Ut will yield a poor approximation of the true 
geodesic distance because it is highly sensitive to errors in mag- 
nitude (see Figures [2] and [6]). The heat method asks for something 
different: it asks only that the gradient Vut points in the right direc- 
tion, i.e., parallel to V(/). Magnitude can safely be ignored since we 
know (from the eikonal equation) that the gradient of the true dis- 
tance function has unit length. We therefore compute the normal- 
ized gradient field X = — Vi^/| Vi^| and find the closest scalar po- 
tential by minimizing | V(/) — Xp, or equivalently, b y solving 
the corresponding Euler-Lagrange equations Ac/) = V-X |Schwarz| 
[T995J . The overall procedure is depicted in Figure [5] 

3.1 Time Discretization 

We discretize the heat equation from step I of Algorithm [T] in time 
using a single backward Euler step for some fixed time t. In prac- 
tice, this means we simply solve the linear equation 



(id — tA)Ut — Uq 



(3) 



over the entire domain M, where id is the identity (here we still 
consider a smooth manifold; spatial discretization is discussed in 




Fig. 6. Left: Varadhan's formula. Right: the heat method. Even for very 
small values of t, simply applying Varadhan's formula does not provide an 
accurate approximation of geodesic distance {top left); for large values of 
t spacing becomes even more uneven {bottom left). Normalizing the gradi- 
ent results in a more accurate solution, as indicated by evenly spaced iso- 
lines { top right), and is also valuable when constructing a smoothed distance 
function {bottom right). 



Section [3^ . Note that backward Euler provides a ma ximum prin- 
ciple, preventing spurious oscillations in the solution j jWade et aT] 
2005 1 . We can get a better understanding of solutions to Eq. ([3J by 
considering the elliptic boundary value problem 



(id - tA)vt = on M\j 
Vt = I on 7 . 



(4) 



which for a point source yields a solution Vt equal to Ut up to a 
multiplicative constant. As established by Varadhan in his proof of 
Eq. ([TJ, Vt also has a close relationship with distance, namely 



lim - 



Vt 
2 



log Vt 



(5) 



This relationship ensures the validity of steps II and III since the 
transformation applied to Vt preserves the direction of the gradient. 

3.2 Spatial Discretization 

In principle the heat method can be applied to any domain with a 
discrete gradient (V), divergence (V-) and Laplace operator (A). 
Here we investigate several possible discretizations. 

3.2. 1 Simplicial Meshes. Let u G M'^' specify a piecewise lin- 
ear function on a triangulated surface. A standard discretization of 
the Laplacian at a vertex i is given by 



(L^ 



1 



(cot aij + cot Pij ) (uj 



where Ai is one third the area of all trian 
gles incident on vertex i, the sum is taken over 
all neighboring vertices j, and aij, /3ij are the 
an gles oppo sing the corresponding edge |Mac 
Neal 1949| . We can express this operation via 
a matrix L A'^Lc, where A e RI^I^I^I is 
a diagonal matrix containing the vertex areas and G 
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Fig. 7. Since the heat method is based on well-established discrete opera- 
tors like the Laplacian, it is easy to adapt to a variety of geometric domains. 
Above: distance on a hippo composed of high-degree nonplanar (and some- 
times nonconvex) polygonal faces. 

is the cotan operator representing the remaining sum. Heat flow 
can then be computed by solving the symmetric positive-definite 
system 



{A - tLc)u - uo 

where uq is a Kronecker delta over 7 (i.e., one 
for source vertices; zero otherwise). The gradient in 
a given triangle can be expressed succinctly as 




2A 



where Af is the area of the face, N is its unit 
normal, is the ith edge vector (oriented counter- ^V/^i^ 
clockwise), and Ui is the value of u at the opposing //x^ { 
vertex. The integrated divergence associated with ver- 
tex i can be written as ^1 

V • X = ^ ^cot6'i(ei • Xj) + cot6'2(e2 • Xj) 



where the sum is taken over incident triangles 2 each with a vec- 
tor Xj, ei and 62 are the two edge vectors of triangle j containing i, 
and Oi , O2 are the opposing angles. If we let d G M'^' be the vector 
of (integrated) divergences of the normalized vector field X, then 
the final distance function is computed by solving the symmetric 
Poisson problem 

Lc0 d. 

Conveniently, this discretization easily generalizes to higher di- 
mensions (e.g., tetrahedral meshes) using well-established discrete 
operators; see for instance [Desbrun et al. 2008) . 

3.2.2 Polygonal Surfaces. For a mesh with (not necessarily 
planar) polygonal faces, we use the polygonal Laplacian defined 
by Alexa and Wardetzky |201 1 1. The only difference in this setting 
is that the gradient of the heat kernel is expressed as a discrete 1- 
form associated with half edges, hence we cannot directly evaluate 
the magnitude of the gradient |Vi^| needed for the normalization 
step (Algorithm [1] step II). To resolve this issue we assume that 
Vu is constant over each face, implying that 



J M 



\Vu\'dA= \Vu\'Af 




where Uf is the vector of heat values in face f , Af is the magnitude 
of the area vector, and Lf is the local (weak) Laplacian. We can 



Fig. 8. The heat method can be applied directly to point clouds that lack 
connectivity information. Left: face scan with holes and noise. Right: kitten 
surface with connectivity removed. Yellow points are close to the source; 
disconnected clusters (in the sense of Liu et al.) receive a constant red value. 

therefore approximate the magnitude of the gradient as 

= ^ujLfUf/Af 

which is used to normalize the 1-form values in the corresponding 
face. The integrated divergence is given by d^Ma where a is the 
normalized gradient, d is the co boundary operator and M is the 
mass matrix for 1 -forms (see [ Al exa and Wardetzky 201 f| for de- 
tails). These operators are applied in steps I-III as usual. Figure [7] 
demonstrates distance computed on an irregular polygonal mesh. 

3.2.3 Point Clouds. For a discrete collection of point samples 
P C of M with no connectivity information, we solve the 
heat equation (step I) using the symmetric point cloud Laplacian 
recently introduced by Liu et al. |2012|, which extends previous 
work of Belkin et al. 1 2009a). In this formulation, the Laplacian 
is represented by Ay^Lpc, where Ay is a diagonal matrix of ap- 
proximate Voronoi areas associated with each p oint, and Lp g is a 
symmetric positive semidefinite matrix (see [Liu et al. 2012) , Sec- 
tion 3.4, for details). 

To compute the vector field X = — Vi^/|Vi^| (step II), we rep- 
resent the function : P ^ R as a height function over ap- 
proximate tangent planes Tp at each point p e P and evaluate 
the gradient of a weighted least squares (WLS) approximation of 
u |Nealen 20 04 1. To compute tangent planes, we use a moving 
least squares (MLS) approximation for simplicity - although other 
choices might be desirable (see Liu et al). The WLS approxima- 
tion of Vu also provides a linear mapping u ^ Du, taking any 
scalar function u to its gradient. To find the best-fit scalar potential 
(step III), we solve the linear, positive-semidefinite Poisson equa- 
tion Lpc(j) — AyX. The distance resulting from this approach 
is depicted in Figure [8] 

Other discretizations are certainly possible (see for instance |Luo| 
et al. 20091); we picked one that was simple to implement in any 
dimension. Note that the computational cost of the heat method de- 
pends primarily on the intrinsic dimension n of M, whereas meth- 
ods based on fast march ing require a grid of the s ame dimension m 
as the ambient space [ Memoli and Sapiro 2001) - this distinction 
is especially important in contexts like machine learning where m 
may be significantly larger than n. 
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3.2.4 Choice of Time Step. The accuracy of the heat method 
reUes in part on the choice of time step t. In the smooth setting, 
Eq. ([5]) suggests that smaller values of t yield better approxima- 
tions of geodesic distance. In the discrete setting we instead dis- 
cover a rather remarkable fact, namely that the limit solution to 
Eq. (|4| is purely a function of the combinatorial distance, inde- 
pendent of how we discretize the Laplacian (see Appendix|A|. The 
main implication of this fact is that - on a fixed mesh - decreasing 
the value of t does not necessarily improve accuracy, even in exact 
arithmetic. (Note, of course, that we can always improve accuracy 
by refining the mesh and decreasing t accordingly.) Moreover, in- 
creasing the value of t past a certain point produces a smoothed ap- 
proximation of geodesic distance (Section [33] l. We therefore seek 
an optimal time step t* that is neither too large nor too small. 

Determining a provably optimal expression for t* is difficult due 
to the grea t complexity of analysis involving the cut locus |NeeT| 
|and Strooc k 2004 1. We instead use a simple estimate that works 
remarkably well in practice, namely t — mh? where h is the mean 
spacing between adjacent nodes and m > is a constant. This 
estimate is motivated by the fact that h'^A is invariant with respect 
to scale and refinement; experiments on a regular grid (Figure [18]) 
suggest that m = 1 is the smallest parameter value that recovers the 
£2 distance, and indeed this value yields near-optimal accuracy for 
a wide variety of irregularly triangulated surfaces, as demonstrated 
in Figure [22] In this paper the time step 



t = h^ 



is therefore used uniformly throughout all tests and examples, ex- 
cept where we explicitly seek a smoothed approximation of dis- 
tance, as in Section [33] 

3.3 Smoothed Distance 

Geodesic distance fails to be smooth at points in the cut locus, i.e., 
points at which there is no unique shortest path to the source - these 
points appear as sharp cusps in the level lines of the distance func- 
tion. Non-smoothness can result in numerical difficulty for applica- 
tions which need to take derivatives of the distance function (e.g., 
level set methods), or may simply be undesirable aesthetically. 

Several distances have been desi gned with smoo thness in mind, 
including diffu sion distance [Coifman and Lafo n 2006 1, commute- 
time distance | Fouss et al. 2007| , and biharmonic distance | |Lipman| 
[et al. 2010] (see the last reference for a more detailed discussion). 
These distances satisfy a number of important properties (smooth- 
ness, isometry-invariance, etc.), but are poor approximations of true 
geodesic distance, as indicated by uneven spacing of isolines (see 
Figure [To| middle). They can also be expensive to evaluate, requir- 
ing either a large number of Laplacian eigenvectors (~ 150 — 200 
in practice) or the solution to a linear system at each vertex. 




Fig. 9. A source on the front of the bunny results in nonsmooth cusps on 
the opposite side. By running heat flow for progressively longer durations 
t, we obtain smoothed approximations of geodesic distance (right). 




Fig. 10. Top row: our smooth approximation of geodesic distance (left) 
and biharmonic distance (middle) both mitigate sharp "cusps" found in 
the exact distance (right), but notice that isoline spacing of the biharmonic 
distance can vary dramatically. Bottom row: biharmonic distance (middle) 
tends to exhibit elliptical level lines near the source, while our smoothed 
distance (left) maintains isotropic circular profiles as seen in the exact dis- 
tance (right). 

In contrast, one can rapidly construct smoothed versions of 
geodesic distance by simply applying the heat method for large val- 
ues of t (Figure [9]). The computational cost remains the same, and 
isolines are evenly spaced for any value of t due to normalization 
(step II). Note that the resulting smooth distance function is isomet- 
rically (but not conformally) invariant since it depends only on the 
intrinsic Laplace-Beltrami operator. 

Interestingly enough, existing smooth distance functions can 
also be understood in terms of time-discrete heat flow. In partic- 
ular, the commute-time distance dc and biharmonic distance 
can be expressed in terms of the harmonic and biharmonic Green's 
functions gc and qb- 

dc{x,yY = gc{x,x) -2gc{x,y)+ gc{y,y), 
ds^x.yY = gB{x,x) -2gB{x,y) + gB{y,y)- 

On a manifold of constant sectional curvature the sum g{x,x) + 
g{y,y) is constant, hence the commute-time and biharmonic dis- 
tances are essentially a scalar multiple of the harmonic and bihar- 
monic Green's functions (respectively), which can be expressed via 
one- and two-step backward Euler approximations of heat flow: 

gc limt^oo(id - tA)^^, 

gB = lim,^oo(id-2tA + t2A2)t^. 

(Here f denotes the pseudoinverse.) Note that for finite t the identity 
operator acts as a regularizer, preventing a logarithmic singularity. 
For spaces with variable curvature, the Green's functions provide 
only an approximation of the corresponding distance functions. 

3.4 Boundary Conditions 

If one is interested in the exact distance, either vanishing Neumann 
or Dirichlet conditions suffice since this choice does not affect the 
behavior of the smooth limit solution (see | von Renesse 2004 1, 
Corollary 2 and [Norris 1997], Theorem 1.1, respectively). Bound- 
ary conditions do however alter the behavior of our smoothed 
geodesic distance (i.e., large t) - Figure [TT] illustrates this behav- 
ior. Although there is no well-defined "correct" behavior for the 
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Fig. 11. Effect of Neumann (top-left), Dirichlet (top-right) and averaged 
(bottom-left) boundary conditions on smoothed distance. Note that aver- 
aged conditions mimic the behavior of the same surface without boundary. 




Fig. 12. For path planning, the behavior of geodesies can be controlled via 
boundary conditions and the integration time t. Top-left: Neumann condi- 
tions encourage boundary adhesion. Top-right: Dirichlet conditions encour- 
age avoidance. Bottom-left: small values of t yield standard straight-line 
geodesies. Bottom-right: large values of t yield more natural trajectories. 

smoothed solution, we advocate the use of averaged boundary con- 
ditions obtained as the mean of the Neumann solution un and the 
Dirichlet solution ud, i-e., u = \{un + ud) - these conditions 
tend to produce isolines that are not substantially influenced by the 
shape of the boundary. The intuition behind this behavior is again 
based on interpreting heat diffusion in terms of random walks: zero 
Dirichlet conditions absorb heat, causing walkers to "fall off" the 
edge of the domain. Neumann conditions prevent heat from flowing 
out of the domain, effectively "reflecting" random walkers. Aver- 
aged boundary conditions mimic the behavior of a domain with- 
out boundary: the number of walkers leaving equals the number of 
walkers returning. Figure[T2]shows how boundary conditions affect 
the behavior of geodesies in a path-planning scenario. 



4. COMPARISON 
4.1 Performance 

A key advantage of the heat 
method is that the linear sys- 
tems in steps (I) and (III) can be 
prefactored. Our implementation 
uses sparse Choleskv f actoriza- 
tion |Chen et al. 2008| , which 
for Pois son- type problems has 
guaranteed sub-quadratic com- 
plexity but in practice scales even 
better [Botsch et al. 2005 1; moreover there is strong evidence 
to suggest that sparse systems arising fr om elliptic PDEs can 
be solved in verv close to linear time [ S chmitz and Ying 2012[ 
Spielman and Teng 2004|. Independent of these issues, the 




amortized cost for problems with a large number of right-hand 
sides is roughly linear, since back substitution can be applied in 
essentially linear time. See inset for a breakdown of relative costs 
in our implementation. 

In terms of absolute performance, a number of factors affect the run 
time of the heat method including the spatial discretization, choice 
of discrete Laplacian, geometric data structures, and so forth. As a 
typi cal exa mple, we compared our simplicial implementation (Sec- 
tion |3.2.1 |) to the first-order fast marching method of Kimmel & 
Sethian |1998 | and the exact algorithm of Mitchell et al 1 19871 
as described by Surazhsky et al |2005 |. In particular we used 
the state-of-the-art fast marching implementation of Peyre and Co- 
hen 12005J and the exact implementation of Kirsanov | Surazhsky] 
|et al. 2005) . The heat method was implemented in ANSI C using 
a simple vertex-face adjacency list. All timings were taken on a 
2.4 GHz Intel Core 2 Duo machine using a single core - Table |l| 
gives timing information. Note that even for a single distance com- 
putation the heat method outperforms fast marching; more impor- 
tantly, updating distance via the heat method for new subsets 7 is 
consistently an order of magnitude faster (or more) than both fast 
marching and the exact algorithm. 

4.2 Accuracy 



We exam ined errors in the heat method, fast marching [Kim mel and 
Sethian 1998], and the exact polyhedral distance | Mitche ll et al. 
1987| , relative to mean edge length h for a variety of triangulated 
surfaces. Figures and [20| illustrate the rate of convergence on 
simple geometries where the smooth geodesic distance can be eas- 
ily obtained. Both fast marching and the heat method appear to ex- 
hibit linear convergence; it is also interesting to note that the exact 
polyhedral distance provides only quadratic convergence. Keeping 
this fact in mind. Table |l| uses the polyhedral distance as a base- 
line for comparison on more complicated geometries - here MAX 
is the maximum error as a percentage of mesh diameter and MiN 
is the mean relative error at each vertex (a convention introduced 
in [Surazhsky et al. 20051 ). Note that fast marching tends to achieve 
a smaller maximum error, whereas the heat method does better on 
average. Figure [T4]gives a visual comparison of accuracy; the only 
notable discrepancy is a slight smoothing at sharp cusps; Figure p3] 
indicates that this phenomenon does not interfere with the extrac- 
tion of the cut locus - here we simply threshold the magnitude of 
A(/). Figure [21] plots the maximum violation of metric properties - 
both the heat method and fast marching exhibit small approxima- 
tion errors that vanish under refinement. Even for the smoothed dis- 
tance (jn » 1) the triangle inequality is violated only for highly 
degenerate geodesic triangles, i.e., all three vertices on a common 
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Fig. 13. Meshes used to test performance and accuracy (see Table[T}. Left to right: Bunny, Isis, Horse, Bimba, Aphrodite, Lion, Ramses. 



Table L Comparison with fast marching and exact polyhedral distance. Best speed/accuracy in bold; speedup in orange. 



Model 


Triangles 




Heat Method 






Fast Marching 


Exact 




Precompute 


Solve 


Max Error 


Mean Error 


Time 


Max Error 


Mean Error 


Time 


Bunny 


28k 


0.21s 


0.01s (28x) 


3.22% 


1.12% 


0.28s 


1.06% 


1.15% 


0.95s 


Isis 


93k 


0.73s 


0.05s (21x) 


1.19% 


0.55% 


1.06s 


0.60% 


0.76% 


5.61s 


Horse 


96k 


0.74s 


0.05s (20x) 


1.18% 


0.42% 


1.00s 


0.74% 


0.66% 


6.42s 


Kitten 


106k 


1.13s 


0.06s (22x) 


0.78% 


0.43% 


1.29s 


0.47% 


0.55% 


11.18s 


Bimba 


149k 


1.79s 


0.09s (29x) 


1.92% 


0.73% 


2.62s 


0.63% 


0.69% 


13.55s 


Aphrodite 


205k 


2.66s 


0.12s (47x) 


1.20% 


0.46% 


5.58s 


0.58% 


0.59% 


25.74s 


Lion 


353k 


5.25s 


0.24s (24x) 


1.92% 


0.84% 


10.92s 


0.68% 


0.67% 


22.33s 


Ramses 


1.6M 


63.4s 


1.45s (68x) 


0.49% 


0.24% 


98.11s 


0.29% 


0.35% 


268.87s 



geodesic. (In contrast, smoothed distances discussed in Section [2] 
satisfy metric properties exactly, but cannot be used to obtain the 
true geometric distance.) Overall, the heat method exhibits errors 
of the same magnitude and rate of convergence as fast marching (at 
lower computational cost) and is likely suitable for any application 
where fast marching is presently used. 

The accuracy of the heat method might be further improved 
by consideri ng alternative s patial discre tizations (see for in- 
stance IBelk in et al. 2009b[ [Hildebrandt and Polthier 201 1|), 
though again one should note that even the exact polyhedral dis- 
tance yields only an 0(/i^) approximation. In the case of the fast 
marching method, accuracy is determined by the choice of update 
rule. A number of highly accurate update rules h ave been devel- 
oped in the case of regular grids {e.g., HJ WENO | Jiang and Peng| 
[1997 1), but fewer options are available on irregular domains such 
as triangle meshes, the predominant choice being the first-order up- 
date rule of Kimmel and Sethian 1 1998 1. Finally, the approximate 
algorithm of Surazhsky et al. provides an interesting comparison 
since it is on par with fast marching in terms of performance and 
produces more accurate results (see [Suraz hsky et al. 2005) , Table 
1). Similar to fast marching, however, it does not take advantage of 
precomputation and therefore exhibits a significantly higher amor- 
tized cost than the heat method; it is also limited to triangle meshes. 




Fig. 14. Visual comparison of accuracy. Left: exact geodesic distance. Us- 
ing default parameters, the heat method (middle) and fast marching (right) 
both produce results of comparable accuracy, here within less than 1% of 
the exact distance - see Table|l|for a more detailed comparison. 
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Fig. 15. Medial axis of the hiragana letter "a" extracted by thresholding 
second derivatives of the distance to the boundary. Left: fast marching. 
Right: heat method. 

4.3 Robustness 

Two factors contribute to the robustness of the heat method, namely 
(1) the use of an unconditionally stable implicit time-integration 
scheme and (2) formulation in terms of purely elliptic PDEs. Fig- 
ure [16] verifies that the heat method continues to work well even on 
meshes that are poorly discretized or corrupted by a large amount 
of noise (here modeled as uniform Gaussian noise applied to the 
vertex coordinates). In this case we use a moderately large value 
of t to investigate the behavior of our smoothed distance; similar 
behavior is observed for small t values. Figure [TTjillustrates the ro- 
bustness of the method on a surface with many small holes as well 
as long sliver triangles. 





Fig. 16. Tests of robustness. Left: our smoothed distance (Section [33) ap- 
pears similar on meshes of different resolution. Right: even for meshes with 
severe noise ( top) we recover a good approximation of the distance function 
on the original surface (bottom, visualized on noise -free mesh). 



Fig. 17. Smoothed geodesic distance on an extremely poor triangulation 
with significant noise - note that small holes are essentially ignored. Also 
note good approximation of distance even along thin slivers in the nose. 

5. CONCLUSION 

The heat method is a simple, general method that can be easily 
incorporated into a broad class of algorithms. However, a great 
deal remains to be explored, including an investigation of alterna- 
tive discretizations. Further improvements on the optimal t value 
also provide an interesting avenue for future work, though typically 
the existing estimate already outperforms fast marching in terms of 
mean error (Table Another obvious question is whether a sim- 
ilar transformation can be applied to a larger class of Hamilton- 
Jacobi equations. Finally, weighted distance computation might be 
achieved by simply rescaling the source data. 
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APPENDIX 

A. A VARADHAN FORMULA FOR GRAPHS 

Lemma 1 . Let G — {V,E) be the graph induced by the spar- 
sity pattern of any real symmetric matrix A, and consider the linear 
system 

(/ - tA)ut = 5 

where L is the identity, 5 is a Kronecker delta at a source vertex 
u ^V, and t > is a real parameter. Then generically 

= lim 

t^o logt 



where G Nq is the graph distance (i.e., number of edges) be- 
tween each vertex v and the source vertex u. 

Proof. Let a be the operator norm of A. Then for t < 1/cr the 
matrix B :— L — tA has an inverse and the solution Ut is given 
by the convergent Neumann series X^^^q i^^^^- Let v e V he 3. 
vertex n edges away from u, and consider the ratio := |s|/|so| 
where Sq :— {t^A'^6)y is the first nonzero term in the sum and 
^ — (X)fcLn+i t^^^^)v is the sum of all remaining terms. Noting 
that \s\ < Er=„+i*'ll^''5|l < ELn+i*''^', we get 

r+ < ^fi—(j _ ^ 

t^{A^6)y 1-ta' 

where the constant c := a^~^-^/{A^6)v does not depend on t. We 
therefore have limt^o = 0, i-^-^ only the first term sq is sig- 
nifcant as t goes to zero. But logso = nlogt + \og{A^6)v is 
dominated by the first term as t goes to zero, hence \og{ut)v/ log t 
approaches the number of edges n. □ 

Numerical experiments such as those depicted in Figure[T8]agree 
with this analysis. 



Fig. 1 8 . Isolines of log ut / log t computed in exact arithmetic on a regular 
grid with unit spacing (h = 1). As predicted by Lemma [T] the solution 
approaches the combinatorial distance as t goes to zero. 
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Fig. 19. convergence of distance functions on the unit sphere with re- 
spect to mean edge length. As a baseline for comparison, we use the exact 
distance function (p{x,y) = cos~-'^(a:: • y). Linear and quadratic conver- 
gence are plotted as dashed lines for reference; note that even the exact 
polyhedral distance converges only quadratically. 




L°° error 



L°° error 



heat method 

fast marching 

exact polyhedral 




Fig. 21. Numerical approximations of geodesic distance exhibit small vi- 
olations of metric properties that vanish under refinement. Here we exam- 
ine errors in symmetry (top left) and the triangle inequality (top right) by 
checking all pairs or triples (respectively) of vertices on the Stanford bunny 
and plotting the worst violation as a percent of mesh diameter. Linear con- 
vergence is plotted as a dashed line for reference. Bottom right: violation 
of triangle inequality occurs only for degenerate geodesic triangles, i.e., all 
three vertices along a common geodesic. Fixing the first two vertices, we 
plot those in violation in red. Bottom left: percent of vertices in violation; 
letting t = mh?, each curve corresponds to a value of m sampled from the 
range [1,100]. 



Fig. 20. Convergence of geodesic distance on the torus at four different 
test points. Error is the absolute value of the difference between the numeri- 
cal value and the exact (smooth) distance; linear and quadratic convergence 
are plotted as dashed lines for reference. Right: test points visualized on the 
torus; dark blue lines are geodesic circles computed via Clairaut's relation. 
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Mean Error 




Fig. 22. Mean percent error as a function of m, where t = mh? . Each 
curve corresponds to a data set from Table [l| Notice that in most examples 
m = 1 (dashed line) is close to the optimal parameter value (blue dots) and 
yields mean error below 1%. 
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